/*====================== MiscIrelandCoordConverter.m ========================*/
/* MiscIrelandCoordConverter class supports the calculations required for
   converting between World and the Irish version of the Universal Transverse
   Mercator coordinates.

   There is only one instance ever, so unless changes are made, this class
   is NON REENTRANT.

   For information on the underlying mathematics, refer to:

	 1- Ordinance Survey Information, "Transverse Mercator Projection,
	    Constants, Formula and Methods", March 1983

	 2- Ordinance Survey, "Tables for the Transverse Mercator Projection
	    of Ireland", 1953, reprinted 1971

NOTES:
	- Irish document uses old tables and some constants to keep
	  values small. This has no effect here.

	- UK document has IIIA and XIIA equation, ie extra correction terms,
	  which are not used on Irish Grid and are left out here.

	- UK document cancels the -'s in V and VI and the -'s in the E
	  equation, which the Irish Grid does not do. This has no effect.

	- Irish document last term in VI is an extra one not present in UK
	  but is included here.

	- There is an erratum to the Irish document that gives the new
	  values of a and b. Old values are in the IrishGridOldUTMConstants,
	  new ones are in the IrishGridUTMConstants object.

	- Accuracy is difficult to compare with the example values in the
	  document because they are derived from tables with interpolation.
	  Our calculated values are within .1 meter of the one from the
	  table, and it is probably due to the fact that the tabular examples
	  do not use "second differences". In other words, the calculated
	  values are probably more correct to more places than the book
	  example values.

IMPORTANT: These equations are accurate to within 1 millimeter. Calculations
	   requiring greater accuracy must use formulae in:

		Redfern, JCB, "Transverse Mercator Formulae", 1948,
		Empire Survey Review, 9(69) pg318-322

	    Extra decimal places are stored only for the purpose of 
	    slowing error propagation that affects the numbers at the
	    millimeter scale, not because they are meaningful in and of 
	    themselves.

   DMA Release 0.8, Copyright @1993 by Genesis Project, Ltd. All Rights
   Reserved. For further information on terms and conditions see:
		Documentation/GISKit/Agreements-Legal-README

HISTORY
12-Mar-93  Dale Amon at GPL
	   Split UTMGrid into constants and convertor parts.
22-Feb-93  Dale Amon at GPL
	   Created.
*/

#import <math.h>
#import <misckit/miscgiskit.h>

@implementation MiscIrelandCoordConverter


/*---------------------------------------------------------------------------*/
/* Calculate latitude and longitude given grid N and E. Accurate to 1 mm.
   Uses constants:
	lambda0 = longitude of grid origin
*/

- (void) utmToWorld
{	double	E,N;
	double	phiPrime;
	double	y,y2,y3,y4,y5,y6;
	double	nu2,nu3,nu4,nu5,nu7;
	double	tanPhiPrime,tan2PhiPrime,tan4PhiPrime,tan6PhiPrime;
	double	secPhiPrime;
	double	VII, VIII,IX, X, XI, XII;

	/* conversion is driven by the source constants */
	xlate = srcConstants;

	E = src[MISC_EASTINGS];
	N = src[MISC_NORTHINGS];

	phiPrime = [self calcPhiPrime: N];

	/* precalculate nu, rho and etaSqrd and powers of nu */
	[self blackboardCalc: phiPrime];
	nu2   = nu*nu;
	nu3   = nu2*nu;
	nu4   = nu2*nu2;
	nu5   = nu3*nu2;
	nu7   = nu4*nu3;

	/* precalculate all the trig values */

	tanPhiPrime  = tan(phiPrime);
	tan2PhiPrime = tanPhiPrime  * tanPhiPrime;
	tan4PhiPrime = tan2PhiPrime * tan2PhiPrime;
	tan6PhiPrime = tan4PhiPrime * tan2PhiPrime;
	secPhiPrime  = 1/cos(phiPrime);

	/* precalculate all the powers of delta E */
	y  = E - xlate->E0;
	y2 = y*y;
	y3 = y2*y;
	y4 = y2*y2;
	y5 = y3*y2;
	y6 = y3*y3;

	/* Calculate the terms used by the latitude and longitude equations */

	VII  = tanPhiPrime/(2.0*rho*nu);
	VIII = tanPhiPrime/(24.0*rho*nu3) *
		(5.0 + 3.0*tan2PhiPrime + etaSqrd - 9.0*etaSqrd*tan2PhiPrime);
	IX   = tanPhiPrime/(720.0*rho*nu5) *
		(61.0 + 90.0*tan2PhiPrime + 45.0*tan4PhiPrime);

	X    = secPhiPrime/nu;
	XI   = secPhiPrime/(6.0*nu3) * (nu/rho + 2.0*tan2PhiPrime);
	XII  = secPhiPrime/(120.0*nu5) *
		(5.0 + 28.0*tan2PhiPrime + 24.0*tan4PhiPrime);

	/* and finally, we give you the latitude and the longitude */	
	dst[MISC_LATITUDE]  = phiPrime     - y2*VII + y4*VIII + y6*IX;
	dst[MISC_LONGITUDE] = xlate->lambda0 + y*X - y3*XI + y5*XII;
	dst[MISC_ALTITUDE]  = src[MISC_ELEVATION];
}


/*---------------------------------------------------------------------------*/
/* Given World Coordinates, calculate the local grid coordinates in
   a UTM projection.
   Uses constants:
	N0,E0 = True origin offset from grid False origin
	lambda0 = longitude of True origin. 
*/

- (void) worldToUTM
{	double	phi,lambda;
	double	cosPhi,cos2Phi,cos3Phi,cos5Phi;
	double	tanPhi,tan2Phi,tan4Phi;
	double	P1, P2, P3, P4, P5;
	double	I,II,III,IV,V,VI;

	/* conversion is driven by the destination constants */
	xlate = dstConstants;

	/* assumes the internal storage format is double precision radians */
	phi = src[MISC_LATITUDE]; lambda = src[MISC_LONGITUDE];

	/* Calculate the Developed Arc of a Meridian from phi to the
	   True Origin. */
	[self calcM: phi];

	/* calculate nu, rho and etaSqrd */
	[self blackboardCalc: phi];

	/* precalculate all the trig values that are not on blackboard */
	cosPhi  = cos(phi);
	cos2Phi = cosPhi*cosPhi;
	cos3Phi = cos2Phi*cosPhi;
	cos5Phi = cos2Phi*cos3Phi;

	tanPhi  = tan(phi);
	tan2Phi = tanPhi  * tanPhi;
	tan4Phi = tan2Phi * tan2Phi;

	/* Precalculate all the powers of delta lambda */
	P1 = lambda - xlate->lambda0;
	P2 = P1 * P1;
	P3 = P2 * P1;
	P4 = P2 * P2;
	P5 = P3 * P2;

	/* Calculate the terms used by the Easting and Northing equations */

	I   = M + xlate->N0;
	II  = nu/2.0   * sinPhi * cosPhi; 
	III = nu/24.0  * sinPhi * cos3Phi * (5.0 - tan2Phi + 9.0 * etaSqrd);
	IV  = nu       * cosPhi;
	V   = nu/6.0   * cos3Phi * (nu/rho - tan2Phi);

	/* The Irish equation has one extra term in VI that the UK
	   and UTM to not have. */
	VI  = nu/120.0 * cos5Phi * 
		(5.0 - 18.0*tan2Phi + tan4Phi + 14.0*etaSqrd - 
		 58.0 * etaSqrd * tan2Phi +  2.0 * etaSqrd * tan4Phi);

	/* And finally, we calculate the local grid coordinates */
	dst[MISC_NORTHINGS] = I         + P2*II + P4*III;
	dst[MISC_EASTINGS]  = xlate->E0 + P1*IV + P3*V   + P5*VI;
	dst[MISC_ELEVATION] = src[MISC_ALTITUDE];
 }


/*===========================================================================*/
/* Class methods */
/*===========================================================================*/
/* Only one converter of this type is every needed. Of course if we got
   into a really big multiprocess agora system there might be a case
   for multiple converters. Since this object is shared by many, it
   doesn't particularly matter what zone it is allocated from.

   Note that we use a secret allocation routine. The parent class blocks the
   normal alloc's because it also creates only a single instance.
*/

static id theIrelandCoordConverter = nil;

+ new
{
    if (!theIrelandCoordConverter)
	theIrelandCoordConverter = [[super superalloc] init];
    return theIrelandCoordConverter;
}


/*===========================================================================*/
/* Initialization methods */
/*===========================================================================*/
/* We add a list of all the services which we are able to provide.
   Note that there is only one MiscIrelandCoordConverter ever created. Once 
   initialized the same object is given to all comers and can not be destroyed.

   Simple copies when the class is the same is okay because we know that
   by definition the UTM constants are be the same.

*/

- init
  { Class	world, utm, oldutm;

    [super init];
    world     = [MiscWorldCoord         class];
    utm       = [MiscIrelandUTMCoord    class];
    oldutm    = [MiscIrelandOldUTMCoord class];

    [self addService: @selector(utmToWorld) convertsFrom: utm    to: world]; 
    [self addService: @selector(worldToUTM) convertsFrom: world  to: utm]; 
    [self addService: @selector(utmToWorld) convertsFrom: oldutm to: world]; 
    [self addService: @selector(worldToUTM) convertsFrom: world  to: oldutm]; 
    [self addService: @selector(copyCoords) convertsFrom: world  to: world]; 
    [self addService:[self fastCopySelector] convertsFrom: utm   to: utm]; 
    [self addService:[self fastCopySelector] convertsFrom:oldutm to: oldutm];

	return self;
  }


/*===========================================================================*/
/* Archive methods */
/*===========================================================================*/

- finishUnarchiving
 {
	[self free];
	return [MiscIrelandCoordConverter new];
 }


@end
